Solutions to Prelim 1 Review Questions
QUESTION 2 ------------------------------------------------------
"M Questions" from Insight. You can check these on your own
using Matlab. Ask a course staff member if you have questions.
QUESTION 3 ------------------------------------------------------
amt = input('Enter the Income: ');
if (amt > 60000)
tax = (25/100)*30000+(30/100)*30000+(40/100)*(amt-60000);
elseif (amt > 30000)
tax = (25/100)*30000+(30/100)*(amt-30000);
else
tax = (25/100)*(amt);
end
fprintf('The tax to be paid for an income of %d is %d',amt,tax);
QUESTION 4 -----------------------------------------------------
% Sec 2.1 Problem 8
% The Euler constant
disp(' n E_n');
disp('-------------------------');
kMax = 100;
step = 100;
% Initialization
total = 0; % Keep track of the summation
for n=1:kMax*step
% Calculate the sum
total = total + 1/n;
if rem(n,step) == 0
% Calculate E_n
E_n = total - log(n);
fprintf(' %5d %.8f\n', n, E_n)
end
end
-----------------------------------------------------
% Sec 2.2 Problem 8
% Two sequences that approximate pi
n = 0;
a = 6/sqrt(3)*(-1)^n/3^n/(2*n+1);
fprintf('a%d = %.7f\n',0, a)
while abs(a - pi) > .000001
n = n + 1;
a = a + 6/sqrt(3)*(-1)^n/3^n/(2*n+1);
fprintf('a%d = %.7f\n',n, a)
end
n = 0;
b = 16 * (-1)^n / 5^(2*n+1) / (2*n+1) - 4 * (-1)^n / 239^(2*n + 1);
fprintf('\nb%d = %.7f\n',0, b)
while abs(b - pi) > .000001
n = n + 1;
b = b + 16 * (-1)^n / 5^(2*n+1) / (2*n+1) - 4 * (-1)^n / 239^(2*n + 1);
fprintf('b%d = %.7f\n',n, b)
end
-----------------------------------------------------
% Section 3.1 Problem 8
% Draw ASCII figure
% *
% * *
% * *
% * *
% * *
% * *
% * *
% * *
% *
n= input('Input n, n>3: ');
% Top half
for r= 1:n
% A general row has spaces, star, spaces, star, \n
for c= 1:n-r
fprintf(' ');
end
fprintf('*');
if (r > 1)
for c = 1:2*(r-1)-1
fprintf(' ');
end
fprintf('*');
end
fprintf('\n')
end
% Bottom half
for r= n-1:-1:1
for c= 1:n-r
fprintf(' ');
end
fprintf('*');
if (r > 1)
for c = 1:2*(r-1)-1
fprintf(' ');
end
fprintf('*');
end
fprintf('\n')
end
-----------------------------------------------------
%Section 3.1 Problem 10
%Print a sequence a(2), a(3), ..., a(n) where n is smallest integer
%such that abs(a(n-1)-a(n))<=.01.
% Initialization
tol= 0.01; % Stopping tolerance
n= 1;
a_current= 0.5; % a(n), i.e., a(1)
a_previous= 0; % a(n-1)
while abs(a_previous - a_current) > tol
% Update
n= n+1;
a_previous= a_current;
% Calculate current term a(n), which is a summation
a_current = 0;
for j = 1:n^2
a_current = a_current + (n/(n^2+j^2));
end
fprintf('%4d %10.6f\n', n, a_current)
end
-----------------------------------------------------
% Section 3.2 Problem 7
% Print t_1 to t_26
disp(' n t_n')
disp('--------------------')
n = 1;
while n<=26
% Initialization
current = 1;
k = n;
% Calculate from inside to outside.
% current is the value of the kth squre root
while k>1 || (n==1 && k>0)
% Update...
current = sqrt(1+k*current);
k = k-1;
end
fprintf(' %2d %10.8f \n',n,current)
% Update...
n = n+1;
end
QUESTION 5 ------------------------------------------------------
% Assume variables nBig and nSm are given.
% Outer loop to iterate from nBig to nSm (while-loop used here;
% you can use for loop).
while nBig >= nSm
% Check nBig to see if it is prime
d=2; % divisor
% Iterate until first proper divisor is found
while (mod(nBig,d) ~= 0);
d = d + 1;
end
if (d == nBig)
fprintf('%d is a prime\n',nBig);
else
fprintf('%d\n',nBig);
end
nBig = nBig � 1;
end
QUESTION 6 -----------------------------------------------------
% Example program for computing the mode of a sequence of integers.
disp('Determine mode of a set of nonnegative integers.')
disp('Use a negative number to quit.')
previous = -1; % Previous number seen
count = 1; % Count of current number
mode = -1; % Mode seen so far
modeCount = 0; % Count for mode so far
number = input('Give me a number: ');
while number >= 0 % Quit when negative number is encountered
if number == previous % same run, so increment count
count = count + 1;
else % new run, so reset count
count = 1;
end
if count > modeCount
mode = number;
modeCount = count;
end
previous = number;
number = input('Another number not smaller than the previous: ');
end
if mode == -1
disp('Mode is undefined')
else
fprintf('Mode is %d which occurred %d times.\n', mode, modeCount)
end
QUESTION 7 -----------------------------------------------------
function seconds = hms2s(h, m, s)
% Convert a time expressed in hours, minutes, seconds to a time in seconds.
% h2ms( h, m, s ) returns 3600*h + 60*m + s
seconds = 3600*h+60*m+s;
%-------------------------------------
% Question 7, cont'd
function [h, m, s] = s2hms( seconds )
% Convert a time in seconds to a time in hrs, minutes, seconds.
% seconds = total number of seconds
% h = maximum number of hours in seconds
% m = maximum number of minutes in seconds, after max hours are removed
% s = number of seconds remaining after max hours and max minutes removed
h = floor( seconds/3600 );
m = floor( rem( seconds, 3600 )/60 );
s = rem( rem(seconds, 3600), 60 );
%-------------------------------------
% Question 7, cont'd
% Script to demonstrate functions hms2s and s2hms
sec = hms2s( 7, 45, 13 )
[h, m, s ] = s2hms( 3682 )
QUESTION 8 -----------------------------------------------------
% Insight Ch4 Sec1 #5
% Solicits input for values a, b,c and prints the length of a:c:b and the
% distance from the last vector component to b.
a = input('a = ');
b = input('b = ');
c = input('c = ');
v = a:c:b;
fprintf('The length of a:c:b is %d.\n', length(v))
fprintf('The distance from v(end) to b is %f.\n', v(end)-b)
% Note: When used as an index, the keyword end is the same as the
% length of the vector. So the distance from the last vector component to
% b can be written as v(length(v)) - b
%-------------------------------------
% Insight Ch4 Sec1 #6
x = 1:10:100;
while length(x) < 1000
x = [x x(1) x];
end
fprintf('length(x) = %d.\n',length(x));
% Note: If the length of x in iteration k is N, then after iteration k+1
% the length of x is 2*N + 1. At the start length(x) = 10 so as the loop
% progresses the length of x is
% 10, 21, 43, 87, 175, 351, 703, 1407.
%-------------------------------------
% Insight Ch4 Sec1 #11
% Generate a length-10 vector according to a recursion formula.
n = input('n = ');
f = zeros(1,10);
for k=1:10
if k==1 % Starting value
f(k) = n;
elseif rem(f(k-1),2)==1 % f(k-1) is odd
f(k) = 3*f(k-1) + 1;
else % f(k-1) is even
f(k) = f(k-1)/2;
end
end
plot(1:10,f,'*') % Plot the result using the star-marker
%-------------------------------------
% Insight Ch5 Sec1 #9
function Snm = Subset(n, m)
% Snm is the number ways a set of n objects can be partitioned into m nonempty subsets.
% n and m are positive integers.
total= 0;
sgn= 1; % sign of the term
for j= m:-1:1 % count down since the mth term has positive sign (sgn is 1)
total= total + sgn * j^n / factorial(m-j) / factorial(j);
sgn= -sgn;
end
Snm= total;
% Note: Due to round-off errors, Snm may not be an integer. To avoid this
% problem, one can calculate the summation of Snm * m! instead, i.e.,
% every term in the summumation is multiplied by m!.
% This means you actually calculate the summation from j=1 to j=m of
% (-1)^m-j * j^n * m! / ((m-j)! * j!)
% ^^^^^^^^^^^^^^^^^^
% This part is the binomial coefficient B(m,j)
% which can be calculated as bchoosek(m,j).
% Refer to problems P5.1.6 and P5.1.7.
%
% Below is the alternative solution:
%
% total = 0;
% sgn = 1; % sign of the term
% for j=m:-1:1
% total = total + sgn * nchoosek(m, j) * j^n;
% sgn = -sgn;
% end
% Snm = total / factorial(m);
%-------------------------------------
% Insight Ch5 Sec1 #9, cont'd
% Print the table of Sigma(n,m)
% print the header line
disp(' n S_n^(1) S_n^(2) S_n^(3) S_n^(4) S_n^(5) S_n^(6)')
disp('---------------------------------------------------------------------')
% print the values
for n=1:10
% string of the line to be printed
line = sprintf('%3d ', n);
for m=1:6
val= Subset(n,m); % val may be non-integer due to round-off error
val= round(val);
line = sprintf('%s%6d ', line, val);
end
disp(line)
end
%-------------------------------------
% Insight Ch6 Sec1 #15
% Simulate selecting two points on the unit circle (uniformly random)
% and calculate the probability that the connecting chord is longer than 1.
N = 10000; % number of trials
num = 0; % number of times that the the connecting chord is longer than 1
for k=1:N
% generate the random points on the unit circle
theta1 = rand(1) * 2*pi;
theta2 = rand(1) * 2*pi;
cx1 = cos(theta1);
cy1 = sin(theta1);
cx2 = cos(theta2);
cy2 = sin(theta2);
% test if the connecting chord is longer than 1
if (cx1 - cx2)^2 + (cy1 - cy2)^2 > 1
num = num + 1;
end
end
fprintf('Probability that connecting chord is longer than 1 is %f\n', num/N)
%-------------------------------------
% Insight Ch6 Sec1 #19a
function p = ProbG(L,R)
% p is the area under the Gaussian function (bell curve) between L and R
% initialize...
N= 10000; % Number of trials
count = 0; % Number of points under the curve
% count the number of random points that lie beneath the curve
for i = 1:N
% Random point in the rectangle bounded by x=L, x=R, y=0, and y=1
x= rand(1)*(R-L)+L;
y= rand(1);
% If the point is under the curve, i.e. f(x), then count as a "hit."
if y <= (2*pi)^(-1/2) * exp(-x^2/2)
count=count+1;
end
end
% Monte Carlo estimate of area under the curve
p= (count/N)*(R-L);
%-------------------------------------
% Insight Ch6 Sec2 #8
function [x, y, z] = RandomWalk3D(n)
% Random walk in 3D. x,y,z are vectors representing the path.
% The kth step has the coordinates (x(k),y(k),z(k)).
%initialize coordinates and counter
xc=0; yc=0; zc=0; k=0;
%while not on the boundary, take a random step
while abs(xc)<n && abs(yc)<n && abs(zc)<n
r=rand(1);
if r <= 1/6
xc = xc + 1;
elseif r <= 2/6
xc = xc - 1;
elseif r <= 3/6
yc = yc + 1;
elseif r <= 4/6
yc = yc - 1;
elseif r <= 5/6
zc = zc + 1;
else
zc = zc - 1;
end
%update position arrays
k=k+1; x(k)=xc; y(k)=yc; z(k)=zc;
end
%-------------------------------------
% Insight Ch6 Sec2 #8, cont'd
% Estimate the average number of hops required for the robot to
% reach the boundary in a 3D random walk.
clc
disp(' n Average #steps to Boundary')
disp('----------------------------------')
nTrials = 100;
for n = 5:5:50
s = 0;
for k=1:nTrials
[x,y,z] = RandomWalk3D(n);
s = s + length(x);
end
ave = s/nTrials;
fprintf(' %3d %8.3f\n',n,ave)
end
fprintf('\n\n(Results based on %d trials)\n',nTrials)